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ABSTRACT 

We propose a formalism for estimating the skewness and angular power spectrum of 
a general Cosmic Microwave Background data set. We use the Edgeworth Expansion 
to define a non-Gaussian likelihood function that takes into account the anisotropic 
nature of the noise and the incompleteness of the sky coverage. The formalism is then 
applied to estimate the skewness of the publicly available 4 year Cosmic Background 
Explorer (COBE) Differential Microwave Radiometer data. We find that the data is 
consistent with a Gaussian skewness, and with isotropy. Inclusion of non Gaussian 
degrees of freedom has essentially no effect on estimates of the power spectrum, if 
each C£ is regarded as a separate parameter or if the angular power spectrum is 
parametrized in terms of an amplitude (Q) and spectral index (n). Fixing the value of 
the angular power spectrum at its maxiumum likelihood estimate, the best fit skewness 
is S" = 6.5 lb 6.0 X 10^(;uK)^; marginalizing over Q the estimate of the skewness is 
S = 6.5 =b 8.4 X 10^(/iK)^ and marginalizing over n one has S = 6.5 ± 8.5 x 10^(/iK)^. 

Subject headings: Cosmology: theory - observation - cosmic microwave background: 
cosmic microwave background 



1. Introduction 

The main assumption in current analysis of Cosmic Microwave Background (CMB) data 
is that the perturbations responsible for structure formation in the universe are Gaussian. 
The theoretical prejudice is that the origin of these fluctuations was a period of Inflation: the 
vacuum fluctuations of the inflaton field on small scales (which satisfy Gaussian statistics) were 
linearly amplified and stretched to super horizon scales through the superluminal expansion of 
the universe. Although there are proposals for non-Gaussian fluctuations within the framework 
of inflation (|Srednicki 1993| , pangui et al 1994 [Peebles 199"7| , [Linde fc Mukhanov 19"9^ ), the 



paradigm is suficiently compelling that the assumption of Gaussianity has been unquestioned 
when constructing methods of data analysis. 

Recently, a number of groups have detected non-Gaussianity in the publicly available 4 
year COBE DMR data. [Ferreira, Magueijo k, Gorski (1998)] found that the COBE normahzed 



bispectrum was inconsistent with that of a Gaussian random field and were unable to atribute 
it to systematic effects or foreground contaminants. Pando, Valls-Gabaud & Fang (1998), using 



a statistic defined in terms of the wavelet transform of the COBE data, found evidence for 
non-Gaussianity localized in the northern hemisphere and Novikov, Feldman and Shandarin 



(1998) identified deviations from Gaussianity using topological measures; Bromley & Tegmark 



have confirmed these results. Banday, Zaroubi &: Gorski (1999) have suggested that the culprit for 



this non-Gaussianity, is a systematic effect in the time ordered data used to construct the publicly 
available COBE maps. 

Possibly the most important result to come out of the analysis of the 4 year COBE data is 
an estimate of the angular power spectrum of the CMB; it can be used to normalize theories of 
large scale structure which when then compared to measures of clustering on other scales can be 
assessed as to their ability to describe our universe. In all estimates of the angular spectrum, 
the fluctuations in the CMB are assumed to be Gaussian. Given that the publicly available 
COBE 4 year data is indeed non-Gaussian, there is a possibility that all estimates of the angular 
power spectrum are biased in some way. In this letter we present a new estimate of the angular 



power spectrum using an approximate likelihood function which includes non-Gaussian degrees of 
freedom. In section we introduce the Edgeworth expansion; in section we apply it to the 4 
year COBE DMR data and in section ^ we discuss our findings. 

2. The Formalism 

In most of what follows we will use the notation of [Gorski (1994)| . A map of the CMB as 
rendered by the COBE DMR instrument, Ap, can be expanded in spherical harmonics 
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where 1^^, are real, orthonormal spherical harmonics, wf are the filter coefficients of the DMR 
beam pattern (Wright et al. 1994), and vJ'f^^ are the filter coefficients used to model smoothing 
due to pixelization; afj^^^ , af^ and af^ are the harmonic coefficients of the CMB anisotropy, 
Galactic emission and instrument noise (from now one we will use the single index notation for 
the spherical harmonic labels, i.e. 1 = 1"^ + (. + 1 + m); iip is the unit vector which points at pixel 
p. A major aim of the DMR analysis is to obtain a good estimate of afjj^, separating out the 
effects due to the Galactic contamination and reciever noise. This can be done by modeling the 



Galactic emission, using information at different frequencies Bennet et al (1992) , and identifing 



which pixels in Ap are contaminated; one then excises these pixels from the map. 

Introducing a Galactic cut in the data renders the Y^^s non-orthonormal. Gorski (1994) 



proposed a solution to this problem: he identified a new basis, ^j(np) which is orthonormal on the 
cut sphere. The aj and the Yi{np can be related to the harmonic coeffecients, Cj and functions in 
the new basis through c = L-^a and V'(%) = L^^Y(np). We have introduced vector notation for 
quantities with index i; L is an triangular matrix so that Cj is a combination of Oj with j > i. It is 
straightforward to see that Ap = c ip{np). 

If the sky is assumed to be Gaussian, then the probability distribution of Fourier amplitudes 
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Here Ccmb = {ccmbCcmb) = ^'^(^cmbslcmb)^^ (^cmbslcmb) = diag(a2) and the noise 
covariance matrix is Cn = ^pix SpG{cutsfcj/} '^p''A(%)V'"^(%) (where Up is the RMS noise in pixel 
p and ^pix is the sohd angle of the pixel). The Bayesian approach is, given the c to find the 
Cg = (a^f£^)) which minimize the Likelihood function, C{C£\C) = G{c\Ci) (where we have 
assumed uniform priors for both c and Ci). 

The Edgeworth expansion permits a systematic extension of the Likelihood function to 
include non-Gaussian degrees of freedom ( [Kendall Sz Stuart 1977 , Chambers 1967, McCullagh 



1984, Juszkiewicz et al 1995, Amendola 1996|); in this letter we shall assume that the non- 



Gaussianity manifests itself as different levels of skewness, 5q,, in sets of pixels, Va- The Edgeworth 
expansion for the skewness is 



F{C^,S^\c) = G{c\C,)[l + Y,\s:£8l + Y.ll^SpS^ll E^pV] 
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where defining W = ip{np){CcMB + Cn) ip (n^') and Ap = Wpp/Api one has: 
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The Edgeworth expansion contains a number of undesirable properties. Seen as a function of 
the data, for fixed parameters, a series truncation is in general not a density (not positive definite) 
and nor unimodal (not a single maximum). These properties follow from the general form of the 
hermite polynomials. For fixed data, as a function of the parameters, any truncation still suffers 
from the same, and other equally serious, problems. We work with a prescription which allows us 
to bypass this inconvenience. We rewrite the likelihood as 



F(x) = G(x) 
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which is equivalent to (^) up to second order. The last factor vanishes should there not be any 
correlations between pixels. We then expand the logarithm of the likelihood in a power series: 

logF = logG + Ef E^.^ + eI E [^^(^P^)1 

+ EE^ E ^^%'-^l^l') (6) 

One can check that exponentiation of (g) differs from (y) in terms higher than second order. 
However exponentiation of (^) leads to a density. In fact the likelihood as a function of s is 
now a Gaussian distribution. Our expansion is formally the same as the Edgeworth expansion; 
however any truncation of a series of the form (o) is very different from a similar truncation of 



the Edgeworth series, and is better behaved. This trick is described in [6.19] of Kendall & Stuart 



1977. One may check that if the data is nearly Gaussian (in the sense that the sample cumulants 



are small), and if there is no noise, then our estimator is the trivial estimator J2p{^p)/^- O^ the 
contrary (^) fails to display a maximum in most circumstances. 

The goal is then to minimize — logF, given by Equations (^) and (^), in terms of Ci and Sa- 
A few comments are in order. Firstly we have derived the most general Edgeworth expansion in 
terms of skewness; a priori it would be natural to atribute one value of the skewness to all the 
pixels in the cut sky (so there is only one a and Va is the set of all pixels in the map). However, 
a more conservative attitude is to relax statistical isotropy and assume that different parts of the 
map will have different levels of skewness. In Section (|3|) we will estimate the skewness of the 
full sky and the skewness of the northern and southern hemispheres. Secondly, the Edgeworth 
expansion is valid in the regime of weak non-Gaussianity. Given that there is evidence of strong 
non-Gaussianity in the COBE one might expect that we are beyond the regime of applicability of 
the formalism. However one must bear in mind that the non-Gaussianity which has been identified 
has been "localized" in harmonic space; statistical tests in pixel space have been consistent with 
the hypothesis the 4 year COBE DMR data is a realization of a Gaussian random field. 



3. Results 

In this section we apply the forniahsm we have developed to the 4 year COBE DMR data. 
We will be using the inverse noise variance weighted, coadded maps of the 53A, 53B, 90A and 
90B COBE-DMR channels at resolution 6, pixelized in the galactic frame. We use the extended 



galactic cut of Banday et al 1997 , and Bennet et al 199(: to remove most of the emission from the 



plane of the Galaxy. We project out the components of the map corresponding to tpi with i = 1,4 
which effectively removes the monopole and dipole in the process and keep 961 coefficients. 

The first questions we are interested in addressing is if the assumption of zero skewness leads 
to a bias in the estimate of the angular power spectrum. We have found that the individual 
estimates of the Ci are negligably affected by introducing skewness; the introducion of the 
Edgeworth expansion introduces a roughly constant renormalization of the likelihood near the 
peak and consequently the Ci estimates and its errorbars are affect by only few percent. It is 
convenient to parametrize the angular power spectrum in terms of a normalization, Q = y -^(i2, 
and a spectral index, n such that 
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In Figure 1 we plot the likelihood contours for Q and n, marginalizing over the skewness. Any 
modification to the Gaussian estimate is sufficiently small that the Gaussian likelihood contours 
and our estimates are indistinguishable. This is a promising result: weakening the assumption of 
Gaussianity does not change in any way existing estimates of the angular power spectrum. 

We now focus on the value and distribution of the skewness. Fixing the various Ci at the 
maximum we find that likelihood as a function of the the skewness takes the form given in Figure 2 
(solid line). By construction it is Gaussian with a peak and l-o" error bars, S = 6.5ib6.0 x 10^ (/xK)^. 
The gaussian point {S = 0) lies well within the estimated probability distibution i.e. our most 
stringent estimate of skewness is consistent with Gaussianity. One can relax the assumption of 
statistical isotropy, and consider the skewness in the South and North hemispheres as separate 
parameters (dashed and dotted lines respectively in Figure 2). We observe that most of the 



statistical significance of our result comes from pixels in the Southern hemisphere. There is no 
obvious reason why this is so: the mean noise is essentially equivalent in the northern and southern 
hemisphere as is its range of values within each hemisphere. 

An important point should be clarified with regards to our estimate: the value we obtain is 
larger larger than what one would expect by applying the naive estimator S = jj — J^n ^n ^^ was 
used in Smoot et al 1994. The latter estimate is only strictly valid if Ap are uncorrelated, not 
the case of the level 6 pixelization of the COBE data. In fact a rough estimate of the number of 
effective (or uncorrelated) pixels indicates that the naive estimate should be a factor of 6-8 times 
smaller than what we get, which is indeed what we find. 

One can attempt to identify correlations between the estimates of Q or n and S by analysing 
the joint likelihhod {Q,S) and {n,S) at the peak of the likelihood or marginalizing over n and 
Q respectively. The results are presented in Figure 3 where we plot the 68% and 95% contour 
levels. Clearly there is very little correlation between the pairs. Marginalizing over n (Figure 3 a) 
or Q ( Figure 3 c) greatly increases the spread of the likelihood in Q and n respecively but has 
essentially the same effect on the effect on the estimated errors of the skewness; In the former 
case the maximum likelihood estimate is S" = 6.5 it 8.4 x 10*^ (/iK)'^ and marginalizing over n one 
has S = 6.5 lb 8.5 x 10^(;uK)^ Marginalizing over both Q and n we obtain or most conservative 
estimate of the skewness, S* = 6.5 it 8.7 x 10*^;uK)^. 

4. Discussion 

We set up a Baysian framework with which to study systematically non Gaussianity. This 
is the natural next step in power spectrum estimation from CMB data sets with more general 
statistics and should be incorporated in future algorithms. We applied this framework to the 
estimate of the skewness of publicly available COBE-DMR 4 year maps. We found that the data 
is consistent with isotropy and Gaussianity as far as the skewness is concerned. In agreement with 
this, the corrections on power spectrum estimation due to introducing the skewness degree of 
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freedom are very small. If we parameterize the power spectrum in terms of all the Ci the effects 
are negligable and one can safely use current Gaussian estimates of the angular power spectrum 
as correct. We then focused on the various possible estimates of the skewness and found our most 
conservative estimate to be S* = —6.5 ± 8.7 x lO'^ fiK^. 

We should perhaps comment on the differences between frequentist and Bayesian studies of 
non Gaussianity. A frequentist approach (such as the one used in Kogut et al (1996)) compares 
the value of a statistic (say skewness) as measured in noisy data with the distribution of what 
should have been measured if the underlying signal were Gaussian. The value of the measured 
statistic is not an estimate for its value for our sky. The Bayesian estimate on the contrary, is an 
estimate of what the signal skewness actually is, given the noise properties of the instrument, and 
what was observed. The Bayesian estimate for the skewness and its frequentist value need not, 
in fact, should not be the same. Nonetheless it is reassuring that both answers agree in that the 
signal's skewness is Gaussian. 

This is the first attempt at estimating the skewness from the publicly available COBE-DMR 4 
year maps; Kogut et al (1996) estimated the pseudocoUapsed and equilateral three point function 
and found these to be also consistent with Gaussianity. Our result reinforces the fact that the 
non-Gaussianity present in COBE maps does not manifest itself in pixel space; the non-Gaussian 
signal is mostly localized in frequency space and only statistics that are sufficiently localized in i 
space will be sensitive to it. The natural next step is to extend the Edgeworth formalism to the 
bispectrum; the complexity of the problem makes it currently numerically intractable for general 
sky coverage and noise covariance matrices. We are currently exploring various regimes where 
such a formalism may be feasible. 
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Fig. 1. — The likelihood as a function of the amplitude Q and spectral index n, marginalized over 
the skewness, 5; the 68%, 95% and 99.7% contours are plotted. The result is indistinguishable 
from the S* = likelihood. 
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Fig. 2. — The Likelihood as a function of skewness, S with the angular power spectrum fixed at its 
maximum likelihood estimates; the dashed line (dotted) line corresponds to using only the northern 
(southern) hemisphere. The solid line corresponds to using both. 
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Fig. 3. — The likelihoods as a function of amplitude, Q and skewness, S, a) marginalized over 
the spectral index, n and b) with n fixed at its maximum likelihood estimate. The likelihoods as 
a function of n and S", c) marginalized over Q and d) with Q fixed at its maximum likelihood 
estimate. 



